Hydra diagnostics

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE,
                      knitr.duplicate.label = "allow")
suppressPackageStartupMessages(library(tidyverse))
library(ggforce)
source("../R/read.report.R")
source("../R/gettables.R")

load data objects

#datafile object
load("../test-data/hydraDataList.rda")

plot data (maybe put at end)

survey_obs <- hydraDataList$observedBiomass %>% 
  mutate(obs = biomass + 1e-07,
         log_obs = log(obs),
         #log_pred = log(pred),
         log_lo = log_obs - 1.96*cv,
         log_hi = log_obs + 1.96*cv,
         obs_lo = exp(log_lo),
         obs_hi = exp(log_hi),
         species = hydraDataList$speciesList[species])

survey indices of abundance

for (surv in unique(survey_obs$survey)) {
  cat("##### Survey #", surv, " \n")
  p1 <- survey_obs %>% 
  filter(survey == surv) %>% 
  ggplot() +
  aes(x= year, y = log_obs, group = species, col=factor(survey)) +
  geom_errorbar(aes(ymin = log_lo, ymax = log_hi)) +
  geom_point() +
#  geom_line(aes(x=year, y=log_pred), col = "blue") +
  facet_wrap(~species, scales = "free_y") +
  theme_bw() +
  guides(species = "None")
  print(p1)
  cat(" \n\n")
  }

Survey # 1

Survey # 2

catch_obs <- hydraDataList$observedCatch %>% 
  mutate(obs = catch + 1e-07,
         log_obs = log(obs),
         #log_pred = log(pred),
         log_lo = log_obs - 1.96*cv,
         log_hi = log_obs + 1.96*cv,
         obs_lo = exp(log_lo),
         obs_hi = exp(log_hi),
         species = hydraDataList$speciesList[species])

Catch time series

for (fleet in unique(catch_obs$fishery)) {
  cat("##### Fleet #", fleet, " \n")
p3 <- catch_obs %>% 
  filter(fishery == fleet,
         catch > 0) %>% 
  ggplot() +
  aes(x= year, y = log_obs, group = species, col= factor(fishery)) +
  geom_errorbar(aes(ymin = log_lo, ymax = log_hi)) +
  geom_point() +
  #geom_line(aes(x=year, y=log_pred), col = "blue") +
  facet_wrap(~species, scales = "free_y") +
  theme_bw() +
  guides(species = "None")
  print(p3)
  cat(" \n\n")
}

Fleet # 1

Size composition data

Diet composition data

Show moar inputs

load run output

repfile <- "../test-data/hydra_sim.rep"
output<-reptoRlist(repfile)

obs_surveyB <- hydraDataList$observedBiomass
obs_catchB <- hydraDataList$observedCatch

biorows <- dim(obs_surveyB)[1]
catchrows <- dim(obs_catchB)[1]
#
indexfits <- gettables(repfile, biorows, catchrows)

show run overview

Model fits

abundance indices

#survey_obs <- hydraDataList$observedBiomass %>% 
survey_obspred <- indexfits[1][[1]] %>%
  mutate(obs = biomass + 1e-07,
         pred = pred_bio + 1e-07,
         log_obs = log(obs),
         log_pred = log(pred),
         log_lo = log_obs - 1.96*cv,
         log_hi = log_obs + 1.96*cv,
         obs_lo = exp(log_lo),
         obs_hi = exp(log_hi),
         species = hydraDataList$speciesList[species])

fits to survey indices of abundance

for (surv in unique(survey_obspred$survey)) {
  cat("##### Survey #", surv, " \n")
  p1 <- survey_obspred %>% 
  filter(survey == surv) %>% 
  ggplot() +
  aes(x= year, y = log_obs, group = species, col=factor(survey)) +
  geom_errorbar(aes(ymin = log_lo, ymax = log_hi)) +
  geom_point() +
  geom_line(aes(x=year, y=log_pred), col = "blue") +
  facet_wrap(~species, scales = "free_y") +
  theme_bw() +
  guides(species = "None")
  print(p1)
  cat(" \n\n")
  }

Survey # 1

Survey # 2

catch

#catch_obs <- hydraDataList$observedCatch %>% 
catch_obspred <- indexfits[2][[1]] %>% 
  mutate(obs = catch + 1e-07,
         pred = pred_catch + 1e-07,
         log_obs = log(obs),
         log_pred = log(pred),
         log_lo = log_obs - 1.96*cv,
         log_hi = log_obs + 1.96*cv,
         obs_lo = exp(log_lo),
         obs_hi = exp(log_hi),
         species = hydraDataList$speciesList[species])

fits to Catch time series

for (fleet in unique(catch_obspred$fishery)) {
  cat("##### Fleet #", fleet, " \n")
p3 <- catch_obspred %>% 
  filter(fishery == fleet,
         catch > 0) %>% 
  ggplot() +
  aes(x= year, y = log_obs, group = species, col= factor(fishery)) +
  geom_errorbar(aes(ymin = log_lo, ymax = log_hi)) +
  geom_point() +
  geom_line(aes(x=year, y=log_pred), col = "blue") +
  facet_wrap(~species, scales = "free_y") +
  theme_bw() +
  guides(species = "None")
  print(p3)
  cat(" \n\n")
}

Fleet # 1

size comps

obs_survey <- hydraDataList$observedSurvSize %>% tibble()

#colnames(obs_survey) <- c('number','year','species','type','InpN','1','2','3','4','5')  #todo - modify so number of bins is based on size of dataframe
obs_survey <- obs_survey %>% pivot_longer(cols=6:ncol(.), names_to = "lenbin") %>%
  filter(value != -999)%>%
  mutate(lenbin = as.integer(str_remove(lenbin, "sizebin")),
         label = rep("survey",nrow(.)),
         species = hydraDataList$speciesList[species])
#obs_survey$label<-rep(("survey"),each=7189)
#dim(obs_survey) #7189

pred_surv<-output$pred_survey_size
obs_survey$pred_surv<-pred_surv
nll_survey<-output$nll_survey_size
obs_survey$nll_surv<-nll_survey
obs_survey$pearson<-((obs_survey$value-obs_survey$pred_surv)/sqrt(obs_survey$pred_surv))

colnames(obs_survey) <- c('number','year','species','type','InpN','lenbin','obs_value','label',
                          'pred_value','nll','pearson')

#obs_catch<-read.table(tsfile, skip=4109, nrows=608, header=F)
obs_catch <- hydraDataList$observedCatchSize %>% tibble()
#colnames(obs_catch) <- c('number','area','year','species','type','InpN','1','2','3','4','5')
obs_catch<-obs_catch %>% pivot_longer(cols=7:ncol(.), names_to = "lenbin") %>%
  mutate(lenbin = as.integer(str_remove(lenbin, "sizebin")),
         label = rep("catch",nrow(.)),
         species = hydraDataList$speciesList[species]) %>%
  filter(value != -999)
#obs_catch$label<-rep(("catch"),each=1855)
#dim(obs_catch) #1855

pred_catch<-output$pred_catch_size
obs_catch$pred_catch<-pred_catch
nll_catch<-output$nll_catch_size
obs_catch$nll_catch<-nll_catch
obs_catch$pearson<-((obs_catch$value-obs_catch$pred_catch)/sqrt(obs_catch$pred_catch))

colnames(obs_catch) <- c('number','area','year','species','type','InpN','lenbin','obs_value','label',
                         'pred_value','nll','pearson')


#diet_catch<-full_join(obs_catch, obs_survey, by = NULL,copy = FALSE)
diet_catch <- bind_rows(obs_catch, obs_survey)

#### READ OBSERVED (.dat) AND PREDICTED (.rep) DIET PROPORTION VALUES ####
# todo as above, make the data file an argument and have the dimensions be determined from the file rather than hard-coded
#obs_diet<-read.table(tsfile, skip=4724, nrows=4721, header=F)
obs_diet <- hydraDataList$observedSurvDiet %>% tibble()

#obs_diet$label<-rep(("diet"),each=4721)
obs_diet<-obs_diet %>% pivot_longer(cols=6:ncol(.), names_to = "prey") %>%
  mutate(#lenbin = as.integer(str_remove(lenbin, "V")),
         species = hydraDataList$speciesList[species],
         label = rep("diet",nrow(.))) %>%
  filter(value >0 )
#obs_diet$label<-rep(("diet"),each=22810)
#dim(obs_diet) #22810

pred_diet<-output$pred_dietprop
obs_diet$pred_diet<-pred_diet
nll_diet<-output$nll_dietprop
obs_diet$nll_diet<-nll_diet
obs_diet$pearson<-((obs_diet$value-obs_diet$pred_diet)/sqrt(obs_diet$pred_diet))
colnames(obs_diet) <- c('number','year','species','lenbin','InpN','prey','obs_value','label', 'pred_value','nll','pearson')

#mydata<-full_join(diet_catch, obs_diet, by = NULL,copy = FALSE)
#mydata <- bind_rows(diet_catch, obs_diet)
mydata <- diet_catch %>%
  mutate(#prey = replace_na(prey, -99),
         area = replace_na(area, 1),
         type = replace_na(type, 0))
data <- mydata
data$residual<-ifelse(data$pearson<0,"negative","positive")
data$res_abs<-abs(data$pearson)
data<-as.data.frame(data)
#head(data)

#names(data)
#str(data)

#select type of data "label= catch, survey or diet"
temp.catch = data[which(data$label == "catch"),]
temp.surv = data[which(data$label == "survey"),]
#temp.diet = data[which(data$label == "diet"),]

data <- obs_diet #data <- read.csv("outputs/survdiet.csv", header = T)
#data<- select(data, -X)
data$residual<-ifelse(data$pearson<0,"negative","positive")
data$res_abs<-abs(data$pearson)
temp.diet = data[which(data$label == "diet"),]

surveys

#### plot 1 length composition plots by species (survey) ####

# modify to go through each survey
# possibly a function that takes data filtered by species and survey as an argument

long_surv <- temp.surv %>%
  pivot_longer(cols = c("pred_value","obs_value"),
               names_to = c("kind","junk"),names_sep = "_") %>%
  select(-junk)

sp<-1
plot_surv <- list()
especies<-unique(long_surv$species)
  for (sp in especies) {
    
    temp_size<-long_surv %>% filter(species == sp & number==1) %>%
      group_by(year) %>%
      summarize(mu_ss=mean(InpN))
    
    plot_surv[[sp]] <- long_surv %>% filter (species==sp& number==1) %>%
      ggplot() +
      aes(x=lenbin, y = value) +
      geom_line(aes(col = kind)) +
      facet_wrap(~year, dir="v") +
      geom_text(data=temp_size, aes(x = 4.5, y = 0.5, label = mu_ss), size=3) +
      theme(legend.position = "bottom") +
      labs(col="") +
      guides(col = guide_legend(nrow = 1))
    
    #ggsave(paste0(plotdir,"complot_surv_",sp,".jpeg"), plot_surv, width = 10, height = 7, dpi = 300)#, width=3000, height=2500, res=250)
    
  }

Size composition of survey 1 by species

for (sp in especies){
  tmp <- plot_surv[[sp]]
  cat("######", sp, " \n")
  print(tmp)
  cat(" \n\n")
  }

Blue_whiting

Capelin

Green_halibut

Haddock

Long_rough_dab

Mackerel

North_atl_cod

Norwegian_ssh

Polar_cod

Redfish

Saithe

#### plot 1 length composition plots by species (survey) ####

sp<-1
plot_surv <- list()
especies<-unique(long_surv$species)
  for (sp in especies) {
    
    temp_size<-long_surv %>% filter(species == sp & number==2) %>%
      group_by(year) %>%
      summarize(mu_ss=mean(InpN))
    
    if(dim(temp_size)[1] > 0) {
    
    plot_surv[[sp]] <- long_surv %>% filter (species==sp& number==2) %>%
      ggplot() +
      aes(x=lenbin, y = value) +
      geom_line(aes(col = kind)) +
      facet_wrap(~year, dir="v", drop = FALSE) +
      geom_text(data=temp_size, aes(x = 4.5, y = 0.5, label = mu_ss), size=3) +
      theme(legend.position = "bottom") +
      labs(col="") +
      guides(col = guide_legend(nrow = 1))
    
    }
    
    if(dim(temp_size)[1] == 0) plot_surv[[sp]] <- NULL
    
    #ggsave(paste0(plotdir,"complot_surv_",sp,".jpeg"), plot_surv, width = 10, height = 7, dpi = 300)#, width=3000, height=2500, res=250)
    
  }

Size composition of survey 2 by species

for (sp in especies){
  tmp <- plot_surv[[sp]]
  cat("######", sp, " \n")
  print(tmp)
  cat(" \n\n")
  }

Blue_whiting

Capelin

Green_halibut

Haddock

Long_rough_dab

Mackerel

NULL

North_atl_cod

Norwegian_ssh

Polar_cod

Redfish

Saithe

catch

long_catch <- temp.catch %>%
  pivot_longer(cols = c("pred_value","obs_value"),
               names_to = c("kind","junk"),names_sep = "_") %>%
  select(-junk)
#sp<-1
#plotdir <- "outputs/figures/comp_plots/by_species/"

especies<-unique(long_catch$species)
plot_catch <- NULL
for (sp in especies) {

  temp_size<-long_catch%>% filter(species == sp) %>%
    group_by(year) %>%
    summarize(mu_ss=mean(InpN))

  plot_catch[[sp]] <- long_catch %>% filter(species==sp) %>%
    ggplot() +
    aes(x=lenbin, y = value) +
    geom_line(aes(col = kind)) +
    facet_wrap(~year, dir="v") +
    geom_text(data=temp_size, aes(x = 4.8, y = 0.5, label = mu_ss), size=3) +
    theme(legend.position = "bottom") +
    labs(col="") +
    guides(col = guide_legend(nrow = 1))

    ##ggsave(paste0(plotdir,"complot_catch_",sp,".jpeg"), plot_catch, width = 10, height = 7, dpi = 300)#, width=3000, height=2500, res=250)

}

Size composition of catch by species

Blue_whiting

Capelin

Green_halibut

Haddock

Mackerel

North_atl_cod

Norwegian_ssh

Redfish

Saithe

Size composition aggregated over time

Catches

#### plot 2 length-comp plots by species aggregated by year (catch and survey) ####

#plotdir <- "outputs/figures/comp_plots/by_year/"

### catch

temporal2 = numeric()
especie = numeric(); especie = sort(unique(temp.catch$species)) # especies
for(e in 1:length(especie)){
  pos1 = numeric(); pos1 = which(temp.catch$species == especie[e])

  temporal1 = numeric()
  year = numeric(); year = sort(unique(temp.catch$year[pos1])) # ano
  for(y in 1:length(year)){
    pos2 = numeric(); pos2 = which(temp.catch$year[pos1] == year[y])

    N = numeric(); N = unique(temp.catch$InpN[pos1][pos2]) # sample size para el estrato especie ano
    n_ejemplares_obs = numeric(); n_ejemplares_obs = round(N * temp.catch$obs_value[pos1][pos2] , 0) # crea vector de proporcion por el sample size
    n_ejemplares_est = numeric(); n_ejemplares_est = round(N * temp.catch$pred_value[pos1][pos2], 0)

    tt1 = numeric()
    tt1 = data.frame(YEAR = year[y], BIN = temp.catch$lenbin[pos1][pos2], N = N, N_EJEM_OBS = n_ejemplares_obs, N_EJEM_EST = n_ejemplares_est)

    temporal1 = rbind(temporal1, tt1) # guarda la proporcion * sample size para la especie seleccionada y ano
  }

  bin1_obs = numeric(); bin1_obs = sum(temporal1$N_EJEM_OBS[which(temporal1$BIN == 1)])
  bin2_obs = numeric(); bin2_obs = sum(temporal1$N_EJEM_OBS[which(temporal1$BIN == 2)])
  bin3_obs = numeric(); bin3_obs = sum(temporal1$N_EJEM_OBS[which(temporal1$BIN == 3)])
  bin4_obs = numeric(); bin4_obs = sum(temporal1$N_EJEM_OBS[which(temporal1$BIN == 4)])
  bin5_obs = numeric(); bin5_obs = sum(temporal1$N_EJEM_OBS[which(temporal1$BIN == 5)])

  total_obs = sum(bin1_obs, bin2_obs, bin3_obs, bin4_obs, bin5_obs)
  prop_new_obs = numeric(); prop_new_obs = c(bin1_obs, bin2_obs, bin3_obs, bin4_obs, bin5_obs) / total_obs

  bin1_pred = numeric(); bin1_pred = sum(temporal1$N_EJEM_EST[which(temporal1$BIN == 1)])
  bin2_pred = numeric(); bin2_pred = sum(temporal1$N_EJEM_EST[which(temporal1$BIN == 2)])
  bin3_pred = numeric(); bin3_pred = sum(temporal1$N_EJEM_EST[which(temporal1$BIN == 3)])
  bin4_pred = numeric(); bin4_pred = sum(temporal1$N_EJEM_EST[which(temporal1$BIN == 4)])
  bin5_pred = numeric(); bin5_pred = sum(temporal1$N_EJEM_EST[which(temporal1$BIN == 5)])

  total_pred = sum(bin1_pred, bin2_pred, bin3_pred, bin4_pred, bin5_pred)
  prop_new_est = numeric(); prop_new_est = c(bin1_pred, bin2_pred, bin3_pred, bin4_pred, bin5_pred) / total_pred

  tt2 = numeric()
  tt2 = data.frame(ESPECIE = especie[e], BIN = seq(1, 5, 1), PROP_OBS = prop_new_obs, PROP_EST = prop_new_est)

  temporal2 = rbind(temporal2, tt2) # almacena por cada especie
}


long_temp2 <- temporal2 %>%
  pivot_longer(cols = c("PROP_OBS","PROP_EST"),
               names_to = c("kind","junk"),names_sep = " ") %>%
  select(-junk)

complot_year<- long_temp2 %>%
  ggplot() +
  aes(x = BIN, y = value) +
  geom_line(aes(col=kind)) +
  theme(title = element_text(angle = 0, hjust = 0.5, size=15, colour="black")) +
  labs(x="length bin", y="proportion value", title="Length composition aggregated by year catch data") +
  facet_wrap(.~ESPECIE) +
  theme(legend.position = "bottom") +
  labs(col="") +
  guides(col = guide_legend(nrow = 1))

#ggsave(paste0(plotdir,"complot_year_catch.jpeg"), complot_year)#, width=3000, height=2500, res=250)
print(complot_year)

Surveys

temporal2 = numeric()
especie = numeric(); especie = sort(unique(temp.surv$species)) # especies
for(e in 1:length(especie)){
  pos1 = numeric(); pos1 = which(temp.surv$species == especie[e])

  temporal1 = numeric()
  year = numeric(); year = sort(unique(temp.surv$year[pos1])) # ano
  for(y in 1:length(year)){
    pos2 = numeric(); pos2 = which(temp.surv$year[pos1] == year[y])

    N = numeric(); N = unique(temp.surv$InpN[pos1][pos2]) # sample size para el estrato especie ano
    n_ejemplares_obs = numeric(); n_ejemplares_obs = round(N * temp.surv$obs_value[pos1][pos2] , 0) # crea vector de proporcion por el sample size
    n_ejemplares_est = numeric(); n_ejemplares_est = round(N * temp.surv$pred_value[pos1][pos2], 0)

    tt1 = numeric()
    tt1 = data.frame(YEAR = year[y], BIN = temp.surv$lenbin[pos1][pos2], N = N, N_EJEM_OBS = n_ejemplares_obs, N_EJEM_EST = n_ejemplares_est)

    temporal1 = rbind(temporal1, tt1) # guarda la proporcion * sample size para la especie seleccionada y ano
  }

  bin1_obs = numeric(); bin1_obs = sum(temporal1$N_EJEM_OBS[which(temporal1$BIN == 1)])
  bin2_obs = numeric(); bin2_obs = sum(temporal1$N_EJEM_OBS[which(temporal1$BIN == 2)])
  bin3_obs = numeric(); bin3_obs = sum(temporal1$N_EJEM_OBS[which(temporal1$BIN == 3)])
  bin4_obs = numeric(); bin4_obs = sum(temporal1$N_EJEM_OBS[which(temporal1$BIN == 4)])
  bin5_obs = numeric(); bin5_obs = sum(temporal1$N_EJEM_OBS[which(temporal1$BIN == 5)])

  total_obs = sum(bin1_obs, bin2_obs, bin3_obs, bin4_obs, bin5_obs)
  prop_new_obs = numeric(); prop_new_obs = c(bin1_obs, bin2_obs, bin3_obs, bin4_obs, bin5_obs) / total_obs

  bin1_pred = numeric(); bin1_pred = sum(temporal1$N_EJEM_EST[which(temporal1$BIN == 1)])
  bin2_pred = numeric(); bin2_pred = sum(temporal1$N_EJEM_EST[which(temporal1$BIN == 2)])
  bin3_pred = numeric(); bin3_pred = sum(temporal1$N_EJEM_EST[which(temporal1$BIN == 3)])
  bin4_pred = numeric(); bin4_pred = sum(temporal1$N_EJEM_EST[which(temporal1$BIN == 4)])
  bin5_pred = numeric(); bin5_pred = sum(temporal1$N_EJEM_EST[which(temporal1$BIN == 5)])

  total_pred = sum(bin1_pred, bin2_pred, bin3_pred, bin4_pred, bin5_pred)
  prop_new_est = numeric(); prop_new_est = c(bin1_pred, bin2_pred, bin3_pred, bin4_pred, bin5_pred) / total_pred

  tt2 = numeric()
  tt2 = data.frame(ESPECIE = especie[e], BIN = seq(1, 5, 1), PROP_OBS = prop_new_obs, PROP_EST = prop_new_est)

  temporal2 = rbind(temporal2, tt2) # almacena por cada especie
}

long_temp2 <- temporal2 %>%
  pivot_longer(cols = c("PROP_OBS","PROP_EST"),
               names_to = c("kind","junk"),names_sep = " ") %>%
  select(-junk)

complot_year<- long_temp2 %>%
  ggplot() +
  aes(x = BIN, y = value) +
  geom_line(aes(col=kind)) +
  theme(title = element_text(angle = 0, hjust = 0.5, size=15, colour="black")) +
  labs(x="length bin", y="proportion value", title="Length composition aggregated by year survey data") +
  facet_wrap(.~ESPECIE) +
  theme(legend.position = "bottom") +
  labs(col="") +
  guides(col = guide_legend(nrow = 1))

#ggsave(paste0(plotdir,"complot_year_survey.jpeg"), complot_year)#, width=3000, height=2500, res=250)
print(complot_year)

Residuals for size composition data

Catch size composition pearson residuals

#### plot 3 pearson residuals bubble plot ####
### catch, species= 1 to 11
#tiff(paste0(plotdir,"bubbleplot_catch.jpeg"), width=3000, height=2500, res=250)
ggplot(temp.catch, aes(x=year, y=lenbin, size = res_abs, color=factor(residual))) +
  geom_point(alpha=0.7) + theme(title = element_text(angle = 0, hjust = 0.5, size=15, colour="black")) +
  facet_wrap(.~species) + labs(x="year", y="length bin", title="Pearson residuals")

#dev.off()

Survey size composition pearson residuals

### survey, species = 1 to 11
#tiff(paste0(plotdir,"bubbleplot_survey.jpeg"), width=3000, height=2500, res=250)
ggplot(temp.surv, aes(x=year, y=lenbin, size = res_abs, color=factor(residual))) +
  geom_point(alpha=0.7) + theme(title = element_text(angle = 0, hjust = 0.5, size=15, colour="black")) +
  facet_wrap(.~species) + labs(x="year", y="length bin", title="Pearson residuals")

dev.off()

## null device 
##           1

Effective sample size plots

Diet composition data

#plotdir <- "outputs/figures/diet/"

predicted<- as.data.frame(cbind(temp.diet$number, temp.diet$year, temp.diet$species,
                                temp.diet$lenbin, temp.diet$pred_value, temp.diet$prey))
colnames(predicted) <- c('number','year','species','lenbin', 'prop', 'prey')
predicted$type2<-"e" #rep(("e"),each=22810)
predicted$sizefit<- paste0(predicted$lenbin,".",predicted$type2)

observed<- as.data.frame(cbind(temp.diet$number, temp.diet$year, temp.diet$species,
                               temp.diet$lenbin, temp.diet$obs_value, temp.diet$prey))
colnames(observed) <- c('number','year','species','lenbin', 'prop', 'prey')
observed$type2<-"o" #rep(("o"),each=22810)
observed$sizefit<- paste0(observed$lenbin,".",observed$type2)

pred_obs <- bind_rows(predicted, observed)# %>%
  #mutate(sizefit = paste0(lenbin,".",type2))

# similar to above for survey length comps, adjust to do for all surveys
# possibly a function that takes filtered data by species, survey as an argument
### species = 1

sp<-1
nsize <- hydraDataList$Nsizebins
stringbit <- paste0(rep(1:nsize, each=2),c(".o",".e"))
limits_use <- rep("",3*nsize)
breaks_use <- rep(NA,3*nsize)
for (i in 1:nsize) {
  lo <- i*3-2
  hi <- i*3-1
  limits_use[lo:hi] <- paste0(rep(i, 2),c(".o",".e"))
  breaks_use[lo:hi] <- limits_use[lo:hi]
}
especies<-unique(pred_obs$species)
plot_diet <- NULL
for (sp in especies) {

  plot_diet[[sp]] <-  pred_obs %>%
    tibble() %>% 
    filter(species == sp & number == 1) %>%
    mutate(prop = as.numeric(prop)) %>% 
    ggplot(aes(x = sizefit, y = prop, group = type2, fill = prey)) +
    geom_col(position = "fill") +
    scale_x_discrete(limits = limits_use,
                     # c("1.o","1.e","",
                     #            "2.o","2.e","",
                     #            "3.o","3.e","",
                     #            "4.o","4.e","",
                     #            "5.o","5.e",""),
                     breaks = breaks_use, #c("1.o","1.e",NA,
                                # "2.o","2.e",NA,
                                # "3.o","3.e",NA,
                                # "4.o","4.e",NA,
                                # "5.o","5.e",NA),
                     labels = limits_use) + #c("1.o","1.e","",
                                # "2.o","2.e","",
                                # "3.o","3.e","",
                                # "4.o","4.e","",
                                # "5.o","5.e","")) +
    coord_flip() +
    facet_wrap(~as.numeric(year)) +
    theme_bw() +
    labs(x = "size & source (o=observed, e=expected)",
         fill = "prey",
         y = "proportion in diet") +
    scale_fill_brewer(type = "qual", palette = 3)
  # facet_wrap_paginate(~ , ncol = 4, nrow = 5, page = 4)

  #ggsave(paste0(plotdir,"diet_plot_",sp,".jpeg"), plot_diet, width = 10, height = 9, dpi = 300)#, width=3000, height=2500, res=250)

}

Proportions by weight for stomach by predator size bin survey 1

Blue_whiting

Green_halibut

Haddock

Long_rough_dab

North_atl_cod

Norwegian_ssh

Redfish

Saithe

sp<-1
nsize <- hydraDataList$Nsizebins
stringbit <- paste0(rep(1:nsize, each=2),c(".o",".e"))
limits_use <- rep("",3*nsize)
breaks_use <- rep(NA,3*nsize)
for (i in 1:nsize) {
  lo <- i*3-2
  hi <- i*3-1
  limits_use[lo:hi] <- paste0(rep(i, 2),c(".o",".e"))
  breaks_use[lo:hi] <- limits_use[lo:hi]
}
especies<-unique(pred_obs$species)
plot_diet <- NULL
for (sp in especies) {
  
  if(dim(pred_obs %>%
    tibble() %>% 
    filter(species == sp & number == 2))[1] > 0) {

  plot_diet[[sp]] <-  pred_obs %>%
    tibble() %>% 
    filter(species == sp & number == 2) %>%
    mutate(prop = as.numeric(prop)) %>% 
    ggplot(aes(x = sizefit, y = prop, group = type2, fill = prey)) +
    geom_col(position = "fill") +
    scale_x_discrete(limits = limits_use,
                     # c("1.o","1.e","",
                     #            "2.o","2.e","",
                     #            "3.o","3.e","",
                     #            "4.o","4.e","",
                     #            "5.o","5.e",""),
                     breaks = breaks_use, #c("1.o","1.e",NA,
                                # "2.o","2.e",NA,
                                # "3.o","3.e",NA,
                                # "4.o","4.e",NA,
                                # "5.o","5.e",NA),
                     labels = limits_use) + #c("1.o","1.e","",
                                # "2.o","2.e","",
                                # "3.o","3.e","",
                                # "4.o","4.e","",
                                # "5.o","5.e","")) +
    coord_flip() +
    facet_wrap(~as.numeric(year)) +
    theme_bw() +
    labs(x = "size & source (o=observed, e=expected)",
         fill = "prey",
         y = "proportion in diet") +
    scale_fill_brewer(type = "qual", palette = 3)
  # facet_wrap_paginate(~ , ncol = 4, nrow = 5, page = 4)
  }
  
  if(dim(pred_obs %>%
    tibble() %>% 
    filter(species == sp & number == 2))[1] == 0) plot_diet[[sp]] <- NULL

  #ggsave(paste0(plotdir,"diet_plot_",sp,".jpeg"), plot_diet, width = 10, height = 9, dpi = 300)#, width=3000, height=2500, res=250)

}

Proportions by weight for stomach by predator size bin survey 2

Blue_whiting

Green_halibut

Haddock

Long_rough_dab

North_atl_cod

Norwegian_ssh

Redfish

Saithe

Effective sample size plots

Model output (derived quantities)

Growth

Maturity

Feeding stuff

Recruitment

Biomass time series

Fishing mortality time series

Predation mortality time series